
#include "qe_math.h"


#define PI				3.14159265358979323846
#define E  				2.7182818284590452354
#define ln_2 			0.69314718055994530942
#define ln_10 			2.30258509299404568402
#define fabs(a)			((a)>0?(a):(-(a)))
#define maxf(a,b) 		((a) > (b) ? (a) : (b))
#define minf(a,b) 		((a) < (b) ? (a) : (b))

static double F1(double x)
{
    return 1/x;
}

static double F2(double x)
{
    return 1/qe_sqrt(1-x*x);
}

static double eee(double x)
{
    if (x > 1e-3) {
        double ee = eee(x/2);
        return ee*ee;
    }
    else {
		 return 1 + x + x*x/2 + qe_pow(x, 3)/6 + 
		 	qe_pow(x, 4)/24 + qe_pow(x, 5)/120;
	}
}

static double simpson(double a, double b, int flag)
{
    double c = a + (b-a)/2;
    if (flag==1)
        return (F1(a)+4*F1(c)+F1(b))*(b-a)/6;
    else
        return (F2(a)+4*F2(c)+F2(b))*(b-a)/6;
}

static double _asr(double a, double b, double eps, double A, int flag)
{
    double c = a + (b-a)/2;
    double L = simpson(a, c,flag), R = simpson(c, b,flag);
    if(fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15.0;
    return _asr(a, c, eps/2, L, flag) + _asr(c, b, eps/2, R, flag);
}

static double asr(double a, double b, double eps, int flag)
{
    return _asr(a, b, eps, simpson(a, b,flag), flag);
}

double qe_sqrt(double x)
{
	if (x > 100) 
		return 10.0*qe_sqrt(x/100);
    double t = x/8 + 0.5 + 2*x/(4+x);
    int c = 10;
    while(c--) {
        t = (t+x/t)/2;
    }
    return t;
}

double qe_pow(double a, int n)
{
	if (n < 0) 
		return 1/qe_pow(a, -n);
	double res = 1.0;
	while (n) {
		if (n & 1) 
			res *= a;
		a *= a;
		n >>= 1;
	}
	return res;
}

double qe_exp(double x)
{
    if (x < 0)
		return 1/qe_exp(-x);
    int n = (int)x;
    x -= n;
    double e1 = qe_pow(E, n);
    double e2 = eee(x);
    return e1*e2;
}

double qe_ln(double x)
{
    return asr(1, x, 1e-8, 1);
}

double qe_powf(double a, double x)
{
    return qe_exp(x*qe_ln(a));
}

